/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */
package proteoformer;

import com.globalmentor.collections.CharSequenceSuffixTree.CharSequenceNode;

import com.globalmentor.collections.CharSequenceSuffixTree;


import java.io.*;
import java.io.IOException;
import java.util.*;
import uk.ac.ebi.pride.tools.mgf_parser.MgfFile;
import uk.ac.ebi.pride.tools.mgf_parser.model.Ms2Query;
import uk.ac.ebi.pride.tools.pride_spectra_clustering.util.ClusteringSpectrum;
import java.lang.Math.*;
import java.text.*;
import com.globalmentor.collections.CharSequenceSuffixTrees;
//import com.globalmentor.collections.SuffixTree;
//import com.globalmentor.collections.SuffixTrees;

/**
 *
 * @author jonesar
 */
public class TagGenerator {
    
    
    //Scan_3866
    //String mgfFile = "C:/Work/ProteomicsSoftware/TopDown/ProteoFormer/svn/ProteoFormer/resources/Scan_2279.mgf";
    static String mgfFile = "C:/Work/ProteomicsSoftware/TopDown/ProteoFormer/svn/ProteoFormer/resources/bsa_2spectra.mgf";  //For XL test
    //String mgfFile = "C:/Work/ProteomicsSoftware/TopDown/ProteoFormer/svn/ProteoFormer/resources/Scan_5265.mgf";
    //String mgfFile = "C:/Work/ProteomicsSoftware/TopDown/ProteoFormer/svn/ProteoFormer/resources/SBEP_STM_004_02222012_Aragon_msdeconv_edit.mgf";
    //String mgfFile = "C:/Work/ProteomicsSoftware/TopDown/ProteoFormer/svn/ProteoFormer/resources/SBEP_STM_004_02222012_Aragon_msdeconv_edit_small.mgf";
    static  String tagFile = mgfFile.substring(0, mgfFile.lastIndexOf('.')) + "_tags.csv";
    static String resFile = mgfFile.substring(0, mgfFile.lastIndexOf('.')) + "_results.csv";
    static String resFullFile = mgfFile.substring(0, mgfFile.lastIndexOf('.')) + "_results_full.csv";
    FileWriter writerRes;
    FileWriter writerFullRes;
    FileWriter writerTags;
    

    //String tagFile = "C:/Work/ProteomicsSoftware/TopDown/ProteoFormer/svn/ProteoFormer/resources/Scan_2279_out.csv";
    
    
    //String fastaFile = "C:/Work/ProteomicsSoftware/TopDown/ProteoFormer/svn/ProteoFormer/resources/ID_002166_F86E3B2F_random.fasta";
    //String fastaFile = "C:/Work/ProteomicsSoftware/TopDown/ProteoFormer/svn/ProteoFormer/resources/OneProt.fasta";
    //String fastaFile = "C:/Work/ProteomicsSoftware/TopDown/ProteoFormer/svn/ProteoFormer/resources/FourProteins.fasta";
    static String fastaFile = "C:/Work/ProteomicsSoftware/TopDown/ProteoFormer/svn/ProteoFormer/resources/bsa_sequence.fasta"; //For XL test
    //String fastaFile = "C:/Work/ProteomicsSoftware/TopDown/ProteoFormer/svn/ProteoFormer/resources/TwoProteins.fasta";
    String protIndexFile = fastaFile.substring(0, fastaFile.lastIndexOf('.')) + "_index.csv";
    //String protIndexFile = "C:/Work/ProteomicsSoftware/TopDown/ProteoFormer/svn/ProteoFormer/resources/ID_002166_F86E3B2F_random_windows1-5_index.csv";
    String prefixIndexFile = protIndexFile.substring(0, protIndexFile.lastIndexOf('.')) + "_prefix.csv";
    String suffixIndexFile = protIndexFile.substring(0, protIndexFile.lastIndexOf('.')) + "_suffix.csv";
    
    String aaMapTable = "C:/Work/ProteomicsSoftware/TopDown/ProteoFormer/svn/ProteoFormer/resources/aaTables.csv";
    

    private static Double HMASS = 1.007825; //TODO - Note, remember to include PROTON mass for calculating charge
    private static Double PROTONMASS = 1.00727647;
    private static Double OMASS = 15.994915;
    private static Double CTERM = 17.002740;
    private static Double NTERM = 1.007825;
    private HashMap<String, Double> aaMap = new HashMap();
    private HashMap<Double,String> aaMapRev = new HashMap();
    private ArrayList<Double> aaMasses = new ArrayList();
    
    static Boolean doProteins = true;
    private HashMap<String, String[]> accToPeptides = new HashMap();
    private HashMap<String, ArrayList<String>> peptideToAccs = new HashMap();
    private String enzymeRegex = "(?<=[KR])(?!P)";      //Trypsin default
    private int missedCleavages = 2;
    
    private HashMap<Double,String> regexAAMapRev = new HashMap();
    private ArrayList<Double> regexAAMasses = new ArrayList();
    
      
    private HashMap<String,ArrayList<Object[]>> proteinIndex = new HashMap();           //The complete protein tag index
    private HashMap<Double,ArrayList<Object[]>> prefixIndex = new HashMap();            //The index accessed by prefix mass - not yet implemented
    private HashMap<Double,ArrayList<Object[]>> suffixIndex = new HashMap();            //The index accessed by suffix mass
    
    //Algorithm notes
    //Where we build regular expression by examining the spectrum - we then need to decide how this can be used to search the tag index
    
    private HashMap<String,ArrayList<Object[]>> idToSpectrumTags = new HashMap();
    
    // **** Following variables required for processing fasta files ***
    private HashMap<String, String> accToSeq = new HashMap();
    private HashMap<String, String> accToDefline = new HashMap(); 
    //private HashMap<Integer,String> intIDToAccessionMap = new HashMap();
    private ArrayList<String> allProteinIDs = new ArrayList();
    private String accessionRegex = " ";
    //*****************************************
    
    private double parts = 50;
    private double million = 1000000;
    private  double ppmError = parts /million;  //initialising as follows doesn't work? ppmError = 3/100000
    private double tagMatchError = 0.03;
    private double MIN_SCORE_THRESHOLD = 5.0;
    private static int MIN_TAG_LENGTH_FOR_SEARCH = 2;

    
        /**
     * @param args the command line arguments
     */
    public static void main(String[] args) {
        
        
        
        TagGenerator tg = new TagGenerator();
        if(args.length!= 5){
            System.out.println("Required arguments for tag generation: fastaFileLocation mgfFileLocation filePrefixLoc true|false(true=proteins) MIN_TAG_LENGTH");            
        }
        else{
            fastaFile = args[0];
            mgfFile = args[1];
            String prefix = args[2];
            doProteins = Boolean.parseBoolean(args[3]);
            MIN_TAG_LENGTH_FOR_SEARCH = Integer.parseInt(args[4]);
            tagFile =  prefix + "_tags.csv";
            resFile =  prefix + "_results.csv";
            resFullFile =  prefix + "_results_full.csv";            
            
        }
        tg.init(fastaFile,mgfFile,tagFile,resFile,resFullFile,doProteins);
        
    }
    
    private void init(String fasta,String mgf, String tagOut, String resOut, String resFullOut, Boolean isProteinTags){
        //testSuffixTree("MPEPTIDERMLLKUGQ");
        //System.exit(0);
        fillAAMap();
        readFasta(fastaFile);
        try{
            /*
            Refinements needed:
                1. What to do in cases where no good tags are generated - not even single amino acids?
              *     -- Look for shifts corresponding to two amino acids
              *     -- Take a look at some of these spectra and see what's going on?
              *     -- Check if different mass error is needed
              *     -- 
                  
            */
            writerTags = new FileWriter(tagFile);
            writerTags.write("spectrumid,start,end,tag,suffix mass,total mass\n");
            loadSpectraFromFile();
            writerTags.close();
            
            //buildProteinIndexFromDisk(protIndexFile);
            
            if(isProteinTags){
                System.out.print("About to build protein index...");
                buildProteinIndex();
            }
            else{
                System.out.println("Building peptide database, assuming " + enzymeRegex + " and missed=" + missedCleavages);
                buildPeptideIndex();
            }
            System.out.print("...done\n");
        
            
            writerRes = new FileWriter(resFile);
            writerFullRes = new FileWriter(resFullFile);
            
            writerFullRes.write("SpectrumID,ProtAcc,peptide,TAG,Prefix,PrefixDelta,Suffix,SuffixDelta\n");                
            writerRes.write("Spectrum,ProtAcc,Seq,Score,MostCommonPrefixMass,CountPrefixMass,MostCommonSuffixMass,CountSuffixMass\n");

            for(String specID : idToSpectrumTags.keySet()){
                ArrayList<Object[]> spectrumTags = idToSpectrumTags.get(specID);
                System.out.println("Querying spectrum: " + specID + " size:" + spectrumTags.size());

                matchPRMToProteins(spectrumTags,specID);
            }
                        
                writerRes.close();
                writerFullRes.close();
         }
        catch(Exception e){
            e.printStackTrace();
        }
    }
    
    
      
    
    
    private void testSuffixTree(String protSeq){
        CharSequenceSuffixTrees trees = new CharSequenceSuffixTrees();
        trees.printTreeNodes(protSeq);
    }

    
    /*
    private void iterateThroughTree(CharSequenceNode node){
        
        System.out.println("\t" + node.toString());
        for(CharSequenceNode childNode : node.)){
            
        }
    }
    */
    
     private void loadSpectraFromFile() {

        try{
            File myFile = new File(mgfFile);
            MgfFile mgfFile = new MgfFile(myFile);
            Iterator<Ms2Query> it = mgfFile.getMs2QueryIterator();
            Set<String> processedIds = new HashSet<String>();

            while(it.hasNext()) {
                Ms2Query query = it.next();

                // make sure every spectrum is only used once
                if (processedIds.contains(query.getTitle()))
                    continue;

                processedIds.add(query.getTitle());

                // set the intensity to 1 in case it's missing
                if (query.getPeptideIntensity() == null)
                    query.setPeptideIntensity(1.0);


                // throw an exception in case it's missing
                if (query.getPrecursorCharge() == null) {
                    throw new Exception("Spectrum is missing precursor charge.");
                }



                // change the id to title
            // ClusteringSpectrum spectrum = new ClusteringSpectrum(query.getTitle(),query.getPrecursorMZ(),query.getPrecursorIntensity(),query.getPrecursorCharge(),query.getPeakList(),null,query.getMsLevel());

                double dummyPrecursorIntensity = 1.0;
                int msLevel =2;

                if(query.getPeakList() != null){
                    int z = query.getPrecursorCharge();
                    double pepMass = (query.getPrecursorMZ()*z)  - (z * PROTONMASS);
                    ArrayList<Double> prmSpectrum = createPRMSpectrum(query);
                    
                    ArrayList<Object[]> spectrumTags = generateTags(query.getTitle(),prmSpectrum,pepMass);
                    System.out.println("Generated tags for: " + query.getTitle()+ " size:" + spectrumTags.size());
                    idToSpectrumTags.put(query.getTitle(), spectrumTags);
                }
                else{
                    System.out.println("\n\n\tNull peak list: "+ query.getTitle());
                }
            }
        }
        catch(Exception e){
            e.printStackTrace();
        }
        
     }
     
     /*
      * Following MS-Align+ method we create the PRM spectrum
      * Currently set to work at 2dp - we might want to use something more complicated to check for same peak at 2dp, then use higher resolution for most intense peak
      */
     private ArrayList<Double> createPRMSpectrum(Ms2Query query){
         
         
         double pepMZ = query.getPrecursorMZ();
         int z = query.getPrecursorCharge();
         double pepMass = (pepMZ  *z) - (z * PROTONMASS);
         String title = query.getTitle();
         
         //System.out.println("Spectrum " + title + " mz: " + pepMZ + " mass: " + pepMass);
         
         // + CTERM + HMASS;
         
         ArrayList<Double> newSpectrum = new ArrayList();
        
         
         HashMap<Double,Double> massToIntensity4dp = new HashMap();     //This will form the spectrum
         HashMap<Double,Double> massToIntensity2dp = new HashMap();     //This is for checking if we have duplication - use the mass at 4dp with higher intensity, add intensity
         
         Map<Double,Double> peakList= query.getPeakList();
         DecimalFormat df4dp = new DecimalFormat("#.####");
         DecimalFormat df2dp = new DecimalFormat("#.##");
         
         //massToIntensity4dp.put(Double.parseDouble(df4dp.format(0.0)), 50.0); //origin 0.0
         //massToIntensity4dp.put(Double.parseDouble(df4dp.format(CTERM +  HMASS)), 50.0); //origin C term group needed for finding y ion series
         //massToIntensity4dp.put(Double.parseDouble(df4dp.format(pepMass)), 50.0); //origin C term group needed for finding y ion series
         
         massToIntensity2dp.put(Double.parseDouble(df2dp.format(0.0)), 50.0); //origin 0.0
         massToIntensity2dp.put(Double.parseDouble(df2dp.format(CTERM +  HMASS)), 50.0); //origin C term group needed for finding y ion series
         massToIntensity2dp.put(Double.parseDouble(df2dp.format(pepMass)), 50.0); //origin C term group needed for finding y ion series
         
         // newSpectrum.add(0.0); //Add zero i.e. to help see if we can find b1 or y1
         //newSpectrum.add(CTERM +  HMASS); //N-terminal group, for chains of y ions
         //newSpectrum.add(pepMass - CTERM -  HMASS );    //PRM of whole protein  TODO - should this be HMASS or PROTONMASS?
         //newSpectrum.add(pepMass );
         
         //first test for duplicate peaks, and include higher intensity ones
         for(Double mz : peakList.keySet()){
             
             //Double formatMass4dp = Double.parseDouble(df4dp.format(mz-PROTONMASS));
             Double formatMass2dp = Double.parseDouble(df2dp.format(mz-PROTONMASS)); 
             double intensity = Double.parseDouble(df4dp.format(peakList.get(mz)));
             
             if(massToIntensity2dp.containsKey(formatMass2dp)){
                intensity += massToIntensity2dp.get(formatMass2dp); 
             }
             massToIntensity2dp.put(formatMass2dp, intensity);
         }
         
         //Now get the inverse mzs and check we are not duplicating
         for(Double mz: peakList.keySet()){
             
             Double inverseMass = Double.parseDouble(df2dp.format(pepMass - (mz-PROTONMASS)));             
             double intensity = Double.parseDouble(df2dp.format(peakList.get(mz)));
             
             /*
             if(massToIntensity4dp.containsKey(inverseMass)){
                 double prevIntensity = massToIntensity4dp.get(inverseMass);                 
                 if(intensity>prevIntensity){
                     massToIntensity4dp.put(inverseMass, intensity);                     
                 }
             }
             else{
                 massToIntensity4dp.put(inverseMass, intensity); 
             } 
             * 
             */
             
             if(massToIntensity2dp.containsKey(inverseMass)){
                intensity += massToIntensity2dp.get(inverseMass); 
             }
             massToIntensity2dp.put(inverseMass, intensity);
         }
         
         for(Double mass : massToIntensity2dp.keySet()){
             //newSpectrum.add(mz);
             newSpectrum.add((mass)); //Change to make this a mass spectrum rather than MH+ spectrum
         }
         
         Collections.sort(newSpectrum);
         
         /* If we want to print out the PRM spectrum
          
         for(Double mass: newSpectrum){
             System.out.println(mass + " " + massToIntensity2dp.get(mass));
         }
*/
         return newSpectrum;
     }
     
     private void fillAAMap() {

        aaMap.put("A", 71.037114);
        aaMap.put("R", 156.101111);
        aaMap.put("N", 114.042927);
        aaMap.put("D", 115.026943);
        aaMap.put("C", 103.009185);
        aaMap.put("E", 129.042593);
        aaMap.put("Q", 128.058578);
        aaMap.put("G", 57.021464);
        aaMap.put("H", 137.058912);
        aaMap.put("I", 113.084064);
        aaMap.put("L", 113.084064);
        aaMap.put("K", 128.094963);
        aaMap.put("M", 131.040485);
        aaMap.put("F", 147.068414);
        aaMap.put("P", 97.052764);
        aaMap.put("S", 87.032028);
        aaMap.put("T", 101.047679);
        aaMap.put("U", 150.95363);
        aaMap.put("W", 186.079313);
        aaMap.put("Y", 163.06332);
        aaMap.put("V", 99.068414);
        //aaMap.put("[", 0.0);             //N-terminus can be inserted for mods
       // aaMap.put("]", 0.0);  
        
        aaMap.put("J", 113.084064);     //ambiguity code for I or L
        
        /*
        aaMapRev.put(71.037114,"A");
        aaMapRev.put(56.101111,"R");
        aaMapRev.put(114.042927,"N");
        aaMapRev.put(115.026943,"D");
        aaMapRev.put(103.009185,"C");
        aaMapRev.put(129.042593,"E");
        aaMapRev.put(128.058578,"Q");
        aaMapRev.put(57.021464,"G");
        aaMapRev.put(137.058912,"H");
        aaMapRev.put(113.084064,"I");
        aaMapRev.put(113.084064,"L");
        aaMapRev.put(128.094963,"K");
        aaMapRev.put(131.040485,"M");
        aaMapRev.put(147.068414,"F" );
        aaMapRev.put(97.052764,"P");
        aaMapRev.put(87.032028,"S");
        aaMapRev.put(101.047679,"T");
        aaMapRev.put(150.95363,"U" );
        aaMapRev.put(186.079313,"W");
        aaMapRev.put(163.06332,"Y");
        aaMapRev.put(99.068414,"V");
        * 
        */
        
        for(String aa : aaMap.keySet()){
            double mz = aaMap.get(aa);
            aaMapRev.put(mz, aa);
        }
        
        Object[] temp = aaMapRev.keySet().toArray();
        for(Object mz : temp){
            aaMasses.add((Double)mz);
        }
        
        DecimalFormat df = new DecimalFormat("#.####");
        for(int i=0;i<aaMasses.size();i++){
                double firstMass = aaMasses.get(i);
                String firstRes = aaMapRev.get(firstMass);
                
                if(regexAAMapRev.containsKey(firstMass)){                    
                    String newPaired=regexAAMapRev.get(firstMass) + "|" + firstRes;
                    regexAAMapRev.put(firstMass, newPaired);
                }
                else{
                    regexAAMapRev.put(firstMass, firstRes);
                }
                
                for(int j=0;j<aaMasses.size();j++){
                    double secondMass = aaMasses.get(j);
                    String secondRes = aaMapRev.get(secondMass);

                    double pairedMass = Double.parseDouble(df.format(firstMass+secondMass));                
                    String paired=firstRes+secondRes;
                    if(regexAAMapRev.containsKey(pairedMass)){                    
                        String newPaired=regexAAMapRev.get(pairedMass) + "|" + paired;
                        regexAAMapRev.put(pairedMass, newPaired);
                    }
                    else{
                        regexAAMapRev.put(pairedMass, paired);
                    }
                                        
                    for(int k=0;k<aaMasses.size();k++){
                        double thirdMass = aaMasses.get(k);
                        String thirdRes = aaMapRev.get(thirdMass);
                        String triplet = firstRes+secondRes+thirdRes;
                        
                        double tripletMass = Double.parseDouble(df.format(firstMass+secondMass+thirdMass));                
                        if(regexAAMapRev.containsKey(tripletMass)){                    
                            String newRegex=regexAAMapRev.get(tripletMass) + "|" + triplet;
                            regexAAMapRev.put(tripletMass, newRegex);
                        }
                        else{
                            regexAAMapRev.put(tripletMass, triplet);
                        }
                    }                    
                }
            }
        
        /*
        for(int i=0;i<aaMasses.size();i++){
            double outerMass = aaMasses.get(i);
            String outerRes = aaMapRev.get(outerMass);
            
            
            
            for(int j=0;j<aaMasses.size();j++){
                double innerMass = aaMasses.get(j);
                String innerRes = aaMapRev.get(innerMass);
                    //String paired = "["+outerRes  + innerRes + "|" + innerRes+outerRes+"]";
                double pairedMass = Double.parseDouble(df.format(outerMass+innerMass));                
                String paired=outerRes+innerRes;
                if(pairedAAMapRev.containsKey(pairedMass)){                    
                    String newPaired=pairedAAMapRev.get(pairedMass) + "|" + paired;
                    pairedAAMapRev.put(pairedMass, newPaired);
                }
                else{
                    pairedAAMapRev.put(pairedMass, paired);
                }
            }
        }

*/
        temp = regexAAMapRev.keySet().toArray();
        for(Object mz : temp){
           regexAAMasses.add((Double)mz);            
        }
        
        try{
            FileWriter writer = new FileWriter(aaMapTable);
            
            /*
            for(Double mz : aaMapRev.keySet()){
                String aa = aaMapRev.get(mz);
                writer.write(df.format(mz)+","+aa+"\n");
            }
             */
            writer.write("mz,regex\n");
            for(Double mz : regexAAMapRev.keySet()){
                String regex = regexAAMapRev.get(mz);
                writer.write(mz+","+regex+"\n");
            }
            
            
            
            
            writer.close();
        }
        catch(IOException e){
            e.printStackTrace();
        }
        
       
        
        Collections.sort(aaMasses);
        Collections.sort(regexAAMasses);
     }
     
     
     private void readFasta(String inputFasta) {

        try {

            BufferedReader in = new BufferedReader(new FileReader(inputFasta));
            String line;

            String currSequence = "";
            String currProtAcc = null;
            String currDefline = null;

            int recordCounter = 0;
            while ((line = in.readLine()) != null) {
                line = line.replaceAll("\n", "");
                line = line.replaceAll("\r", "");

                if (line.contains(">")) {
                    //Insert previous into hash and reset
                    if (recordCounter != 0) {
                        currSequence = currSequence.replaceAll(" ", "");
                        allProteinIDs.add(currProtAcc);

                        String[] peptides = currSequence.split(enzymeRegex);
                        addPeptidesToMap(currProtAcc, peptides);

                        accToSeq.put(currProtAcc, currSequence);
                        //System.out.println("Inserting:" + currProtAcc + "_" + currSequence);
                        accToDefline.put(currProtAcc, currDefline);
                        //System.out.println("Inserting2:" + currProtAcc + "_" + currDefline);

                        currSequence = "";
                    }

                    line = line.replaceAll(">", "");

                    int splitPos = line.indexOf(accessionRegex);
                    if (splitPos != -1) {
                        currProtAcc = line.substring(0, splitPos);
                        currDefline = line.substring(splitPos + 1);
                    }
                    else {
                        System.out.println("Regular expression not found for split: " + line);
                        System.exit(0);
                    }

                    recordCounter++;
                }
                else {
                    currSequence += line;
                }
            }
            //handle last
            currSequence = currSequence.replaceAll(" ", "");
            accToSeq.put(currProtAcc, currSequence);
            accToDefline.put(currProtAcc, currDefline);

            String[] peptides = currSequence.split(enzymeRegex);
            addPeptidesToMap(currProtAcc, peptides);
            //System.out.println("Put last:" + currProtAcc + " seq: " + currSequence);
            allProteinIDs.add(currProtAcc);

            in.close();
        }
        catch (IOException e) {
            System.err.println("Error: " + e.getMessage());
        }
    }
    
     private void addPeptidesToMap(String protAcc, String[] peptides) {

        ArrayList<String> tempPeptides = new ArrayList();

        if (missedCleavages == 0) {
            accToPeptides.put(protAcc, peptides);
        }
        else {

            for (int i = 0; i < peptides.length; i++) {
                String outerPep = peptides[i];
                tempPeptides.add(outerPep);

                //TODO - test missed cleaves code!
                String tempPeptide = "";
                for (int j = i + 1; j <= i + missedCleavages; j++) {

                    if (j < peptides.length) {
                        tempPeptide += peptides[j];
                        String mc = outerPep + tempPeptide;
                        tempPeptides.add(mc);
                        //System.out.println("Adding mc:" + mc );
                    }
                }
            }

            String[] newPeptides = new String[tempPeptides.size()];
            for (int i = 0; i < tempPeptides.size(); i++) {
                newPeptides[i] = tempPeptides.get(i);
            }

            accToPeptides.put(protAcc, newPeptides);

        }

    }
     
     
     
     private ArrayList<Object[]> generateTags(String specTitle,ArrayList<Double> spectrum, double pepMass){
        
         ArrayList<Object[]> spectrumTags = new ArrayList();
         try{
            //FileWriter file = new FileWriter(tagFile);

            double biggestMass = aaMasses.get(aaMasses.size()-1)+0.5;
            double smallestMass = aaMasses.get(0)-0.5;
            
            HashMap<Double,String> endPosToTag = new HashMap();
            
            Collections.sort(spectrum);
            
            /*
             * Algorithm
             * 
             * Build chains:
             * 
             * 1. Insert end positions (peak 2) into a hashMap, where value is the current chain
             * 2. If new start position is a current end position, create new chain with paired values
             * 
             */

            DecimalFormat df = new DecimalFormat("#.####"); //format to 4dp
            
            //First pass through the spectrum looks for perfect tags
            for(int i = 0; i<spectrum.size();i++){

                Double outerMass = Double.parseDouble(df.format(spectrum.get(i)));
                //System.out.println("o:" + outerMass);
                for(int j = i+1; j<spectrum.size();j++){
                    
                    Double innerMass = Double.parseDouble(df.format(spectrum.get(j)));
                    Double delta = Double.parseDouble(df.format(innerMass - outerMass));
                    //System.out.println("\ti:" + innerMass);

                    Boolean aaFound = false;
                    if(delta > smallestMass && delta < biggestMass){
                        double topRange = delta + (ppmError * innerMass);
                        double bottomRange = delta - (ppmError * innerMass);
                        //System.out.print("\t\tppm: " + ppmError + " Delta: " + delta + " bottom: " + bottomRange + " top: " + topRange+"\n");

                        
                        for(int k=0;k<aaMasses.size();k++){
                            Double aa = aaMasses.get(k);
                            //System.out.print("\t\t\t" + aa);

                            if(aa > bottomRange && aa < topRange){
                                //System.out.print("Found" + aa);
                                aaFound = true;
                                double error = delta - aa;
                                double prefix = outerMass;
                                double suffix = pepMass - innerMass ;
                                String foundAA =  aaMapRev.get(aa);
                                String resLine = foundAA + "," + outerMass + "," + innerMass + "," + delta + "," + error + "," + suffix;
                                
                                if(endPosToTag.containsKey(outerMass)){
                                    String currTag = endPosToTag.get(outerMass)+foundAA;    //Add to chain
                                    endPosToTag.put(innerMass, currTag);
                                    //System.out.println("Adding to tag, inner: " + innerMass + " outer: " + outerMass + " tag: " +currTag);
                                }
                                else{
                                    double prefixMass = outerMass;
                                    /*
                                    if(doProteins){
                                        prefixMass =  outerMass- CTERM -  HMASS; //I don't see why this works but it seems to TODO TO DO
                                    }
                                    else{
                                        prefixMass =  outerMass;
                                    }
*/
                                    endPosToTag.put(innerMass,""+outerMass+","+prefixMass+","+foundAA);
                                    //System.out.println("Found tag, inner: " + innerMass + " outer: " + outerMass + " tag: " +foundAA);
                                }
                                
                                //writerTags.write(resLine +"\n");
                                k=aaMasses.size();
                            }
                        }                                              
                        //System.out.print("\n");
                    }
                }                
            }
            
            //Trace through = prefixMass is outerMZ (i.e. first peak - 18Daltons) - unclear why this was needed
            //Suffix Mass = peptide mass - peak being the end of the tag

            
           // writerTags.write("\n\n");
            for(Double endPosition : endPosToTag.keySet()){
                String tempTag = endPosToTag.get(endPosition);
                
                String temp[] = tempTag.split(",");
                String startPos = temp[0];
                String tag = temp[2];
                double suffixMass = pepMass - endPosition;
                
                writerTags.write(specTitle + ","+startPos+","+endPosition+","+tag+","+suffixMass+","+pepMass+"\n");
                double prefixMass = Double.parseDouble(temp[1]);
                
                //TODO - There is redunancy in the tags which needs to be removed
                if(tag.length()>=MIN_TAG_LENGTH_FOR_SEARCH){
                    Object[] triplet = {tag,prefixMass,suffixMass};
                    spectrumTags.add(triplet);
                }                
            }
            
            //Second pass through the data to generate longer strings, including ambiguity
            
             //Reset if allowed paired amino acids
            biggestMass = regexAAMasses.get(regexAAMasses.size()-1)+0.5;
            smallestMass = regexAAMasses.get(0)-0.5;
            endPosToTag = new HashMap();    //start again
            for(int i = 0; i<spectrum.size();i++){

                Double outerMass = Double.parseDouble(df.format(spectrum.get(i)));
                for(int j = i+1; j<spectrum.size();j++){
                    Double innerMass = Double.parseDouble(df.format(spectrum.get(j)));
                    Double delta = Double.parseDouble(df.format(innerMass - outerMass));


                    Boolean aaFound = false;
                    if(delta > smallestMass && delta < biggestMass){
                        double topRange = delta + (ppmError * innerMass);
                        double bottomRange = delta - (ppmError * innerMass);
                        //System.out.print("ppm: " + ppmError + " Delta: " + delta + " bottom: " + bottomRange + " top: " + topRange);

                        
                        for(int k=0;k<aaMasses.size();k++){
                            Double aa = aaMasses.get(k);
                            //System.out.print("\t" + aa);

                            if(aa > bottomRange && aa < topRange){
                                aaFound = true;
                                double error = delta - aa;
                                double prefix = outerMass;
                                double suffix = pepMass - innerMass;
                                String foundAA =  aaMapRev.get(aa);
                                //String resLine = foundAA + "," + outerMZ + "," + innerMZ + "," + delta + "," + error + "," + suffix;
                                
                                if(endPosToTag.containsKey(outerMass)){
                                    String currTag = endPosToTag.get(outerMass)+foundAA;    //Add to chain
                                    endPosToTag.put(innerMass, currTag);
                                }
                                else{
                                    double prefixMass =  outerMass;
                                    
                                    /*
                                    if(doProteins){
                                        
                                        prefixMass = outerMass- CTERM -  HMASS; //I don't see why this works but it seems to TODO TO DO
                                    }
                                    else{
                                        
                                        prefixMass = outerMass;
                                    }
                                    * 
                                    */
                                    endPosToTag.put(innerMass,""+outerMass+","+prefixMass+","+foundAA);
                                }
                                
                                //writerTags.write(resLine +"\n");
                                k=aaMasses.size();
                            }
                        }
                        
                        /*
                         * Now check for pairs of amino acids if diff between peaks is not a single AA
                         * 
                         */
                        if(!aaFound){
                            for(int k=0;k<regexAAMasses.size();k++){
                                Double aa = regexAAMasses.get(k);
                                //System.out.print("\t" + aa);

                                if(aa > bottomRange && aa < topRange){
                                    double error = delta - aa;
                                    double prefix = outerMass;
                                    double suffix = pepMass - innerMass;
                                    String foundRegex = "["+regexAAMapRev.get(aa)+"]";
                                    //System.out.println(foundPair + "\t" + outerMZ + "\t" + innerMZ+ "\t" + delta + "\t" + error + "\t" + suffix);
                                    k=regexAAMasses.size();
                                    if(endPosToTag.containsKey(outerMass)){
                                        String currTag = endPosToTag.get(outerMass)+foundRegex;    //Add to chain
                                        endPosToTag.put(innerMass, currTag);
                                    }
                                    else{
                                        double prefixMass = outerMass;
                                        /*
                                        if(doProteins){
                                            prefixMass =  outerMass- CTERM -  HMASS; //I don't see why this works but it seems to TODO TO DO
                                        }
                                        else{
                                            prefixMass =  outerMass;
                                        }
                                        */
                                        endPosToTag.put(innerMass,""+outerMass+","+prefixMass+","+foundRegex);
                                    }
                                }
                            }
                        }
                        
                        
                        //System.out.print("\n");
                    }
                }                
            }
            
            for(Double endPosition : endPosToTag.keySet()){
                String tempTag = endPosToTag.get(endPosition);
                String temp[] = tempTag.split(",");
                String startPos = temp[0];
                String tag = "%%"+temp[2];//To denote this is now a regular expression
                double suffixMass = pepMass - endPosition;
                
                writerTags.write(specTitle + ","+startPos+","+endPosition+","+tag+","+suffixMass+","+pepMass+"\n");
                double prefixMass = Double.parseDouble(temp[1]);
                
                //TODO - There is redunancy in the tags which needs to be removed
                if(tag.contains("[")){
                    Object[] triplet = {tag,prefixMass,suffixMass};
                    spectrumTags.add(triplet);
                }                
            }

            /*
            for(int i = 0; i<spectrum.size();i++){

                Double outerMZ = Double.parseDouble(df.format(spectrum.get(i)));
                for(int j = i+1; j<spectrum.size();j++){
                    Double innerMZ = Double.parseDouble(df.format(spectrum.get(j)));
                    Double delta = Double.parseDouble(df.format(innerMZ - outerMZ));

                    if(delta > smallestMass && delta < biggestMass){
                        double topRange = delta + (ppmError * innerMZ);
                        double bottomRange = delta - (ppmError * innerMZ);
                        //System.out.print("ppm: " + ppmError + " Delta: " + delta + " bottom: " + bottomRange + " top: " + topRange);

                        for(int k=0;k<pairedAAMasses.size();k++){
                            Double aa = pairedAAMasses.get(k);
                            //System.out.print("\t" + aa);

                            if(aa > bottomRange && aa < topRange){
                                double error = delta - aa;
                                double prefix = outerMZ;
                                double suffix = pepMass - innerMZ;
                                System.out.println(pairedAAMapRev.get(aa) + "\t" + outerMZ + "\t" + innerMZ+ "\t" + delta + "\t" + error + "\t" + suffix);
                                k=pairedAAMasses.size();
                            }
                        }
                        //System.out.print("\n");
                    }
                }
            }
              */       
            
            
          }
         catch(IOException e){
                e.printStackTrace();
         }
         
         return spectrumTags;
     }
     
     private void buildProteinIndexFromDisk(){
         
           try {

                BufferedReader in = new BufferedReader(new FileReader(protIndexFile));
                String line;
                while ((line = in.readLine()) != null) {
                                            //writer.write(acc+","+tag+","+String.format("%.3g%n",prefixMass)+","+String.format("%.3g%n",suffixMass)+"\n");
                    //Object[] triplet = {acc,prefixMass,suffixMass};
                    
                    String[] temp = line.split(",");
                    String tag = temp[0];
                    
                    ArrayList<Object[]> dataList = new ArrayList();
                    proteinIndex.put(tag, dataList);
                    for(int i=1;i<temp.length;i=i+3){                    
                        String acc = (String)temp[i];
                        double prefixMass = Double.parseDouble(temp[i+1]);
                        double suffixMass = Double.parseDouble(temp[i+2]);
                        Object[] triplet = {acc,prefixMass,suffixMass};
                        dataList.add(triplet);
                    }
                    
                }
                
                in = new BufferedReader(new FileReader(prefixIndexFile));
                
                while ((line = in.readLine()) != null) {
                    String[] temp = line.split(",");
                    Double prefix = Double.parseDouble(temp[0]);
                    //Object[] prefixTriplet = {acc,tag,suffixMass};
                    //Object[] suffixTriplet = {acc,tag,prefixMass};
                    ArrayList<Object[]> dataList = new ArrayList();
                    prefixIndex.put(prefix, dataList);
                    for(int i=1;i<temp.length;i=i+3){                    
                        String acc = (String)temp[i];
                        String tag = (String)temp[i+1];
                        double suffixMass = Double.parseDouble(temp[i+2]);
                        Object[] triplet = {acc,tag,suffixMass};
                        dataList.add(triplet);
                    }                    
                }
                
                in = new BufferedReader(new FileReader(suffixIndexFile));
                
                while ((line = in.readLine()) != null) {
                    String[] temp = line.split(",");
                    Double suffix = Double.parseDouble(temp[0]);
                    //Object[] prefixTriplet = {acc,tag,suffixMass};
                    //Object[] suffixTriplet = {acc,tag,prefixMass};
                    ArrayList<Object[]> dataList = new ArrayList();
                   suffixIndex.put(suffix, dataList);
                    for(int i=1;i<temp.length;i=i+3){                    
                        String acc = (String)temp[i];
                        String tag = (String)temp[i+1];
                        double prefixMass = Double.parseDouble(temp[i+2]);
                        Object[] triplet = {acc,tag,prefixMass};
                        dataList.add(triplet);
                    }                    
                }
                
            }
           catch(Exception e){
               e.printStackTrace();
           }
     } 
     
     private void buildProteinIndex(){
         
         //String protSeq1 = "MAAKIRRDDEVIVLTGKDKGKRGKVKNVLSSGKVIVEGINLVKKHQKPVPALNQPGGIVEKEAAIQVSNVAIFNAATGKADRVGFRFEDGKKVRFFKSNSETIK";
         //String protSeq2 = "AAKIRRDDEVIVLTGKDKGKRGKVKNVLSSGKVIVEGINLVKKHQKPVPALNQPGGIVEKEAAIQVSNVAIFNAATGKADRVGFRFEDGKKVRFFKSNSETIK";
         //String protID2 = "prot2";
         //0.0	            AAKIRRDDEVIVLTGKDKGKRGKVKNVLSSGKVIVEGINLVKKHQKPVPALNQPGGIVEKEAAIQVSNVAIFNAATGKADRVGFRFEDGKKVRFFKSNSE	10547.863012999995
         
         try{
            FileWriter writer = new FileWriter(protIndexFile);
            FileWriter prefixWriter = new FileWriter(prefixIndexFile);
            FileWriter suffixWriter = new FileWriter(suffixIndexFile);
            DecimalFormat df = new DecimalFormat("#.####");
            for(String acc : accToSeq.keySet()){
                String protSeq = accToSeq.get(acc);
                protSeq = protSeq.replaceAll("L", "J");
                protSeq = protSeq.replaceAll("I", "J");   //Replacing all Is and Ls with Js
                //int[] windowSizes = {3,4,5,6,7};
                int[] windowSizes = {1,2,3,4,5};

                for(int windowSize : windowSizes){
                    for(int i=0;i<=protSeq.length()-windowSize;i++){             
                        String prefix = protSeq.substring(0,i);
                        String tag = protSeq.substring(i,i+windowSize);
                        String suffix = protSeq.substring(i+windowSize,protSeq.length());

                        double prefixMass=Double.parseDouble(df.format(getPeptideNeutralMass(prefix)));
                        double suffixMass=Double.parseDouble(df.format(getPeptideNeutralMass(suffix)));

                        Object[] tagTriplet = {acc,prefixMass,suffixMass};
                        Object[] prefixTriplet = {acc,tag,suffixMass};
                        Object[] suffixTriplet = {acc,tag,prefixMass};

                        ArrayList<Object[]> dataList = null;
                        if(!proteinIndex.containsKey(tag)){
                            dataList = new ArrayList();
                        }
                        else{
                            dataList = proteinIndex.get(tag);
                            proteinIndex.put(tag,dataList);
                        }
                        dataList.add(tagTriplet);
                        proteinIndex.put(tag, dataList);
                        
                        //Now the prefix index
                       dataList = null;
                       
                       //To implement this as a key search, the keys will need to go down to 1 or 2dp I guess, will also need to implement some kind of fuzzy matching between tags and sequence
                       //        Also need to remember to implement I/L ambiguity...
                       //        Just try out for one or two difficult spectra to see how it works, then stop for now
                       
                        if(!prefixIndex.containsKey(prefixMass)){
                            dataList = new ArrayList();
                        }
                        else{
                            dataList = prefixIndex.get(prefixMass);
                            prefixIndex.put(prefixMass,dataList);
                        }
                        dataList.add(prefixTriplet);
                        prefixIndex.put(prefixMass, dataList);
                        
                        //Now the suffix index
                        dataList = null;
                        if(!suffixIndex.containsKey(suffixMass)){
                            dataList = new ArrayList();
                        }
                        else{
                            dataList = suffixIndex.get(suffixMass);
                            suffixIndex.put(suffixMass,dataList);
                        }
                        dataList.add(suffixTriplet);
                        suffixIndex.put(suffixMass, dataList);
                        
                    }
                }
            }
            
            
            for(String tag : proteinIndex.keySet()){
                ArrayList<Object[]> allDataForTag = proteinIndex.get(tag);
                //String tagLine = tag + ",";
                writer.write(tag);
                for(Object[] data : allDataForTag){                
                    String acc = (String)data[0];
                    double prefixMass = (double)data[1];
                    double suffixMass = (double)data[2];
                    writer.write(","+acc+","+prefixMass+","+suffixMass);
                } 
                writer.write("\n");
            }
            
            for(Double prefix : prefixIndex.keySet()){
                ArrayList<Object[]> allDataForTag = prefixIndex.get(prefix );
                //String tagLine = tag + ",";
                prefixWriter.write(""+prefix);
                for(Object[] data : allDataForTag){                
                    String acc = (String)data[0];
                    String tag = (String)data[1];
                    double suffixMass = (double)data[2];
                    prefixWriter.write(","+acc+","+tag+","+suffixMass);
                } 
                prefixWriter.write("\n");
            }
            
            for(Double suffix : suffixIndex.keySet()){
                ArrayList<Object[]> allDataForTag = suffixIndex.get(suffix );
                //String tagLine = tag + ",";
                suffixWriter.write(""+suffix);
                for(Object[] data : allDataForTag){                
                    String acc = (String)data[0];
                    String tag = (String)data[1];
                    double prefixMass = (double)data[2];
                    suffixWriter.write(","+acc+","+tag+","+prefixMass);
                } 
                suffixWriter.write("\n");
            }
            
            
            writer.close();
            prefixWriter.close();
            suffixWriter.close();
         }
         catch(Exception e){
            e.printStackTrace();
         }

     }

    
     /*
      * Different option to build the index for individual peptides rather than proteins
      * 
      */
    private void buildPeptideIndex(){
         
         //String protSeq1 = "MAAKIRRDDEVIVLTGKDKGKRGKVKNVLSSGKVIVEGINLVKKHQKPVPALNQPGGIVEKEAAIQVSNVAIFNAATGKADRVGFRFEDGKKVRFFKSNSETIK";
         //String protSeq2 = "AAKIRRDDEVIVLTGKDKGKRGKVKNVLSSGKVIVEGINLVKKHQKPVPALNQPGGIVEKEAAIQVSNVAIFNAATGKADRVGFRFEDGKKVRFFKSNSETIK";
         //String protID2 = "prot2";
         //0.0	            AAKIRRDDEVIVLTGKDKGKRGKVKNVLSSGKVIVEGINLVKKHQKPVPALNQPGGIVEKEAAIQVSNVAIFNAATGKADRVGFRFEDGKKVRFFKSNSE	10547.863012999995
         
         try{
            FileWriter writer = new FileWriter(protIndexFile);
            FileWriter prefixWriter = new FileWriter(prefixIndexFile);
            FileWriter suffixWriter = new FileWriter(suffixIndexFile);
            
            DecimalFormat df = new DecimalFormat("#.####");
            for(String acc : accToSeq.keySet()){
                //String protSeq = accToSeq.get(acc);
                
                for(String peptide : accToPeptides.get(acc)){
                    peptide = peptide.replaceAll("L", "J");
                    peptide = peptide.replaceAll("I", "J");   //Replacing all Is and Ls with Js
                    //int[] windowSizes = {3,4,5,6,7};
                    int[] windowSizes = {1,2,3,4,5};
                    
                    //Store peptide to protein maps
                    ArrayList<String> accsFromThisPep = null;
                    if(peptideToAccs.containsKey(peptide)){
                        accsFromThisPep = peptideToAccs.get(peptide);
                    }
                    else{
                        accsFromThisPep = new ArrayList();
                        peptideToAccs.put(peptide, accsFromThisPep);
                    }
                    accsFromThisPep.add(acc);

                    for(int windowSize : windowSizes){
                        for(int i=0;i<=peptide.length()-windowSize;i++){             
                            String prefix = peptide.substring(0,i);
                            String tag = peptide.substring(i,i+windowSize);
                            String suffix = peptide.substring(i+windowSize,peptide.length());

                            double prefixMass=Double.parseDouble(df.format(getPeptideNeutralMass(prefix)));
                            double suffixMass=Double.parseDouble(df.format(getPeptideNeutralMass(suffix)));

                            Object[] tagTriplet = {peptide,prefixMass,suffixMass};
                            Object[] prefixTriplet = {peptide,tag,suffixMass};
                            Object[] suffixTriplet = {peptide,tag,prefixMass};
        
                            ArrayList<Object[]> dataList = null;
                            if(!proteinIndex.containsKey(tag)){
                                dataList = new ArrayList();
                            }
                            else{
                                dataList = proteinIndex.get(tag);
                                proteinIndex.put(tag,dataList);
                            }
                            dataList.add(tagTriplet);
                            proteinIndex.put(tag, dataList);

                            //Now the prefix index
                        dataList = null;

                        //To implement this as a key search, the keys will need to go down to 1 or 2dp I guess, will also need to implement some kind of fuzzy matching between tags and sequence
                        //        Also need to remember to implement I/L ambiguity...
                        //        Just try out for one or two difficult spectra to see how it works, then stop for now

                            if(!prefixIndex.containsKey(prefixMass)){
                                dataList = new ArrayList();
                            }
                            else{
                                dataList = prefixIndex.get(prefixMass);
                                prefixIndex.put(prefixMass,dataList);
                            }
                            dataList.add(prefixTriplet);
                            prefixIndex.put(prefixMass, dataList);

                            //Now the suffix index
                            dataList = null;
                            if(!suffixIndex.containsKey(suffixMass)){
                                dataList = new ArrayList();
                            }
                            else{
                                dataList = suffixIndex.get(suffixMass);
                                suffixIndex.put(suffixMass,dataList);
                            }
                            dataList.add(suffixTriplet);
                            suffixIndex.put(suffixMass, dataList);

                        }
                    }
                }
            }
            
            
            for(String tag : proteinIndex.keySet()){
                ArrayList<Object[]> allDataForTag = proteinIndex.get(tag);
                //String tagLine = tag + ",";
                writer.write(tag);
                for(Object[] data : allDataForTag){                
                    String peptide = (String)data[0];
                    double prefixMass = (double)data[1];
                    double suffixMass = (double)data[2];
                    writer.write(","+peptide+","+prefixMass+","+suffixMass);
                } 
                writer.write("\n");
            }
            
            for(Double prefix : prefixIndex.keySet()){
                ArrayList<Object[]> allDataForTag = prefixIndex.get(prefix );
                //String tagLine = tag + ",";
                prefixWriter.write(""+prefix);
                for(Object[] data : allDataForTag){                
                    String peptide = (String)data[0];
                    String tag = (String)data[1];
                    double suffixMass = (double)data[2];
                    prefixWriter.write(","+peptide+","+tag+","+suffixMass);
                } 
                prefixWriter.write("\n");
            }
            
            for(Double suffix : suffixIndex.keySet()){
                ArrayList<Object[]> allDataForTag = suffixIndex.get(suffix );
                //String tagLine = tag + ",";
                suffixWriter.write(""+suffix);
                for(Object[] data : allDataForTag){                
                    String peptide = (String)data[0];
                    String tag = (String)data[1];
                    double prefixMass = (double)data[2];
                    suffixWriter.write(","+peptide+","+tag+","+prefixMass);
                } 
                suffixWriter.write("\n");
            }
            
            
            writer.close();
            prefixWriter.close();
            suffixWriter.close();
         }
         catch(IOException e){
            e.printStackTrace();
         }

     }

     
     

        /*
        * Algorithm:
        * 1. Attempt to match in tag 7, 6, ... 3 index
        * 2. If get a tag match:
        *      - Check prefix
        *      - Check suffix
        *      - Collect the delta errors from prefix and suffix match - if same mass is seen more than once, assume this is a systematic shift
        * 3. Build a scoring function, based on tag length, prefix and suffix matches. If same prefix or suffix error match is seen more than once, then assume this is a systematic error and not down weighted
        */
        private void matchPRMToProteins(ArrayList<Object[]> spectrumTags,String specID){

            try{
            
                HashMap<String,ArrayList<Object[]>> mapProtToMatches = new HashMap();

                for(Object[] spectrumTag : spectrumTags){
                    
                    DecimalFormat df = new DecimalFormat("#.###");
                    String tag = (String)spectrumTag[0];
                    tag = tag.replaceAll("L", "J");
                    tag = tag.replaceAll("I", "J");   //Replacing all Is and Ls with Js
                    
                    double prefix = (double)spectrumTag[1];
                    double suffix = (double)spectrumTag[2];

                    //System.out.println("tagsearch:"+tag+"\t"+df.format(prefix)+"\t"+df.format(suffix));

                    if(!tag.contains("%%")){
                    
                        
                        ArrayList<Object[]> tagMatches  = proteinIndex.get(tag);
                        ArrayList<Object[]> revTagMatches = proteinIndex.get(new StringBuilder(tag).reverse().toString());
                        //ArrayList<Object[]> allMatches = new ArrayList();

                        /*
                        if(tagMatches!=null){
                            for(int i=0;i<tagMatches.size();i++){
                                allMatches.add(tagMatches.get(i));
                            }
                        }
                        if(revTagMatches!=null){
                            for(int i=0;i<revTagMatches.size();i++){
                                allMatches.add(revTagMatches.get(i));
                            }
                        }
*/
                        
                        
                        
                        if(tagMatches!=null){
                            for(Object[] matchedTag : tagMatches){
                                //{protID2,prefixMass,suffixMass};
                                String id = (String)matchedTag[0];
                                double prefixDelta = prefix - (double)matchedTag[1];    //These need to be reversed in rev matches
                                double suffixDelta = suffix - (double)matchedTag[2]- CTERM -  HMASS; 
                                
                                if(doProteins){
                                    writerFullRes.write(specID+","+id+",null,"+tag+","+prefix+","+prefixDelta+","+suffix+","+suffixDelta+"\n");
                                }
                                else{
                                    ArrayList<String> protAccs = peptideToAccs.get(id);
                                    String accs = protAccs.toString().replaceAll("\\[", "").replaceAll("\\]", "").replaceAll(",", ";");
                                    writerFullRes.write(specID+","+accs+","+id+","+tag+","+prefix+","+prefixDelta+","+suffix+","+suffixDelta+",Forward\n");
                                }

                                ArrayList<Object[]> dataForProt = null;
                                if(!mapProtToMatches.containsKey(id)){
                                    dataForProt = new ArrayList();
                                    mapProtToMatches.put(id,dataForProt);
                                }
                                else{
                                    dataForProt = mapProtToMatches.get(id);
                                }
                                Object[] tagMatchedWithDeltas = {tag,prefixDelta,suffixDelta};    //TODO at source, can save memory (perhaps?) by moving down to 3 dp?
                                dataForProt.add(tagMatchedWithDeltas);
                            }
                        }
                        if(revTagMatches!=null){
                            for(Object[] matchedTag : revTagMatches){
                                //{protID2,prefixMass,suffixMass};
                                String id = (String)matchedTag[0];
                                double prefixDelta = suffix - (double)matchedTag[1];
                                double suffixDelta = prefix - (double)matchedTag[2]- CTERM -  HMASS;    
                                
                                if(doProteins){
                                    writerFullRes.write(specID+","+id+",null,"+tag+","+prefix+","+prefixDelta+","+suffix+","+suffixDelta+"\n");
                                }
                                else{
                                    ArrayList<String> protAccs = peptideToAccs.get(id);
                                    String accs = protAccs.toString().replaceAll("\\[", "").replaceAll("\\]", "").replaceAll(",", ";");
                                    writerFullRes.write(specID+","+accs+","+id+","+tag+","+prefix+","+prefixDelta+","+suffix+","+suffixDelta+",Rev\n");
                                }

                                ArrayList<Object[]> dataForProt = null;
                                if(!mapProtToMatches.containsKey(id)){
                                    dataForProt = new ArrayList();
                                    mapProtToMatches.put(id,dataForProt);
                                }
                                else{
                                    dataForProt = mapProtToMatches.get(id);
                                }
                                Object[] tagMatchedWithDeltas = {tag,prefixDelta,suffixDelta};    //TODO at source, can save memory (perhaps?) by moving down to 3 dp?
                                dataForProt.add(tagMatchedWithDeltas);
                            }
                        }
                        
                    }
                    else{
                        
                        //System.out.println("Regular expression matching not yet implemented");
                        
                    }

                }

                for(String acc : mapProtToMatches.keySet()){
                    ArrayList<Object[]> tagMatches = mapProtToMatches.get(acc);
                    double protScore = 0.0;

                    String protRes = "";
                    
                    HashMap<Double,Integer> mapPrefixErrorsToCount = new HashMap();
                    HashMap<Double,Integer> mapSuffixErrorsToCount = new HashMap();

                    //First pass to pick up the distributions of errors
                    for(Object[] tagMatch : tagMatches){
                        Boolean oneMassMatched = false;
                        String tag = (String)tagMatch[0];
                        int tagScore = tag.length() * tag.length() * tag.length();
                        double tagErrorScore = 1.0 * tagScore;
                        double prefixMass = (double)tagMatch[1];
                        double suffixMass = (double)tagMatch[2];

                        //TODO The following code is rather clunky in the use of HashMaps and Integers
                        DecimalFormat df = new DecimalFormat("#.#");
                        Double roundedPrefixMass = Double.parseDouble(df.format(prefixMass));
                        Double roundedSuffixMass = Double.parseDouble(df.format(suffixMass));

                        //System.out.println("Errors: " + roundedPrefixMass + " \t" + roundedSuffixMass);
                        
                        Integer prefixErrorCount = null;
                        if(!mapPrefixErrorsToCount.containsKey(roundedPrefixMass)){
                            prefixErrorCount = new Integer(1);                        
                        }
                        else{
                            prefixErrorCount = mapPrefixErrorsToCount.get(roundedPrefixMass);
                            prefixErrorCount = new Integer(prefixErrorCount.intValue() + 1);
                        }
                        mapPrefixErrorsToCount.put(roundedPrefixMass,prefixErrorCount);

                        Integer suffixErrorCount = null;
                        if(!mapSuffixErrorsToCount.containsKey(roundedSuffixMass)){
                            suffixErrorCount = new Integer(1);     
                        }
                        else{
                            suffixErrorCount = mapSuffixErrorsToCount.get(roundedSuffixMass);
                            suffixErrorCount = new Integer(suffixErrorCount.intValue() + 1);
                        }
                        mapSuffixErrorsToCount.put(roundedSuffixMass,suffixErrorCount);

                        if(Math.abs(prefixMass) < tagMatchError && Math.abs(suffixMass) < tagMatchError){
                            tagErrorScore *= 1000;
                        }
                        else if(Math.abs(prefixMass) < tagMatchError || Math.abs(suffixMass) < tagMatchError){
                            tagErrorScore *= 300;
                        }

                        protScore += tagErrorScore;
                    }
                    protScore /= 1000;
                    
                    //Now we work out the most frequent prefix and suffix errors
                    Double mostFreqPrefixError = 0.0;
                    int countMostFreqPrefixError =0;                    
                    for(Double prefixError : mapPrefixErrorsToCount.keySet()){
                        
                        int count = mapPrefixErrorsToCount.get(prefixError).intValue();  
                        
                        if(count>countMostFreqPrefixError){
                            countMostFreqPrefixError = count;
                            mostFreqPrefixError = prefixError;
                        }                        
                    }
                    
                    Double mostFreqSuffixError = 0.0;
                    int countMostFreqSuffixError =0;
                    for(Double suffixError : mapSuffixErrorsToCount.keySet()){
                        int count = mapSuffixErrorsToCount.get(suffixError ).intValue();  
                        if(count>countMostFreqSuffixError){
                            countMostFreqSuffixError = count;
                            mostFreqSuffixError = suffixError ;
                        }                        
                    }
                    String seq = accToSeq.get(acc);
                    if(protScore > MIN_SCORE_THRESHOLD){
                        writerRes.write(specID+","+acc+","+seq + "," + protScore+","+mostFreqPrefixError +","+countMostFreqPrefixError +","+ mostFreqSuffixError+","+countMostFreqSuffixError+"\n");
                    }
                }

            }
            catch(Exception e){
                e.printStackTrace();
            }
            
            
            
        }        
        
    
     
     
    /*
     * Calculates the mass of a peptide based on its amino acids only i.e. a b ion
     *
     */
    private double getPeptideNeutralMass(String peptide) {
        double mass = 0.0;

        peptide = peptide.toUpperCase();
        peptide = peptide.trim();
        //System.out.println("Peptide:" + peptidePlusModString);

        /*
        String[] temp = peptidePlusModString.split("_");
        String peptide = "[" + temp[0] + "]";
        String modString = temp[1];

*/
        
        String[] aas = peptide.split("");
        for (String aa : aas) {
            if (aa.length() == 1) {
                if (aaMap.get(aa) != null) {
                    mass += aaMap.get(aa);
                }
                else {
                    System.out.println("Non-recognized char:" + aa);
                }
            }
        }

        //System.out.println( " mr: " + mr); 

        return mass;

    }


    
}
